function jd=cal2jd(yr,mn,dy)
% CAL2JD  Converts calendar date to Julian date using algorithm
%   from "Practical Ephemeris Calculations" by Oliver Montenbruck
%   (Springer-Verlag, 1989). Uses astronomical year for B.C. dates
%   (2 BC = -1 yr). Non-vectorized version. See also DOY2JD, GPS2JD,
%   JD2CAL, JD2DOW, JD2DOY, JD2GPS, JD2YR, YR2JD.
% Version: 2011-11-13
% Usage:   jd=cal2jd(yr,mn,dy)
% Input:   yr - calendar year (4-digit including century)
%          mn - calendar month
%          dy - calendar day (including factional day)
% Output:  jd - jJulian date

% Copyright (c) 2011, Michael R. Craymer
% All rights reserved.
% Email: mike@craymer.com

if nargin ~= 3
  warning('Incorrect number of input arguments');
  return;
end
if any(mn) < 1 || any(mn) > 12
  warning('Invalid input month');
  return
end
if dy < 1
  if any(mn == 2 & dy > 29) || any(any(mn == [3 5 9 11]) && dy > 30) || any(dy > 31)
    warning('Invalid input day');
    return
  end
end

y = yr;
m = mn;
i = mn <= 2;
y(i) = yr(i) - 1;
m(i) = mn(i) + 12;

date1=4.5+31*(10+12*1582);   % Last day of Julian calendar (1582.10.04 Noon)
date2=15.5+31*(10+12*1582);  % First day of Gregorian calendar (1582.10.15 Noon)
date=dy+31*(mn+12*yr);
b = NaN(size(yr));
b(date <= date1) = -2;
i = date >= date2; 
b(i) = fix(y(i)/400) - fix(y(i)/100);
% if any(isnan(b))
%   warning('cal2jd_vect: Dates between October 5 & 15, 1582 do not exist');
% %   return;
% end

jd = fix(365.25*y) + fix(30.6001*(m+1)) + b + 1720996.5 + dy;
i = y <= 0;
jd(i) = fix(365.25*y(i)-0.75) + fix(30.6001*(m(i)+1)) + b(i) + 1720996.5 + dy(i);
